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Abstract 

We describe a three-dimensional hydrodynamic lattice-gas model of amphophilic fluids. This model 
of the non-equilibrium properties of oil-water-surfactant systems, which is a non-trivial extension of an 
earlier two-dimensional realisation due to Boghosian, Coveney and Emerton El, can be studied effec- 
tively only when it is implemented using high-performance computing and visualisation techniques. We 
describe essential aspects of the model's theoretical basis and computer implementation, and report on 
the phenomenological properties of the model which confirm that it correctly captures binary oil-water 
and surfactant-water behaviour, as well as the complex phase behaviour of ternary amphiphilic fluids. 
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1 Introduction 

Oil and water are immiscible fluids under normal conditions of temperature and pressure. Their phase 
separation behaviour in binary mixtures has been extensively studied both experimentally and theoretically, 
and has become a testbed for many fluid simulation methods |^,|l5[. However, the addition of amphiphile (or 
surfactant) to these systems produces much more complex behaviour. A general review of the equilibrium 
phase behaviour of these fluids and theoretical methods for studying them was provided by Gompper and 
Schick @. The complex phase behaviour of binary and ternary amphiphilic fluids arises as a result of 
the physical and chemical properties of amphiphilic molecules. Amphiphile molecules possess a hydrophylic 
head and a hydrophobic tail. As a result it is usually energetically favourable for the amphiphile molecules 
to be adsorbed at oil-water interfaces, effectively reducing the interfacial tension. Complex morphologies, 
termed mesophases, occur in both binary (surfactant and water or oil) and ternary amphiphilic systems. In 
general these structures depend on the concentration and chemical nature of the surfactant molecules (length 
of hydrocarbon tail, size of head group, and so on) as well as on the temperature. The equilibrium phase 
diagrams for these systems have been obtained from experimental investigation |l6| ] . Of considerable interest 
are the microemulsion phases; an example is the very low surface tension which exist between the "middle" 
microemulsion phase and the bulk oil and water phases fl4|| , but there are many other fascinating phenomena 
including the viscoelastic properties of wormlike micelles and the formation of vesicles |^,[l4|] . Their intrinsic 
complexity and wide application make these systems appropriate for detailed scientific inquiry. Although the 
equilibrium phase behaviour of these systems is well understood, relatively little work has been done on their 
non-equilibrium properties. Much of the work which has addressed these dynamic properties has been based 
on molecular dynamics methods Jl8|, [l9|, p2| . However, such atomistic approaches are too computationally 
demanding to allow investigation of the important large-time dynamics and extended spatial structures of 
these systems. Indeed, a large part of the fascination of amphiphilic fluids is related to the dependence 
of their macroscopic properties on the underlying molecular and mesoscopic dynamics, which calls for the 
development of techniques which can efficiently bridge the length and time scale gaps from micro to macro. 
In the present paper we report on the formulation and implementation of a three-dimensional hydrody- 
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namic lattice gas model, which has been extensively studied in two dimensions P, p|~^2[ ■ Compared with 
molecular dynamics the computational simplicity of lattice gases makes them an ideal method for modelling 
complex fluid behaviour from the mesoscale upwards. They have been used extensively for modelling hydro- 
dynamics since Frisch, Hasslacher, and Pomeau |i3fl , and Wolfram ^lj showed that it is possible to simulate 
the incompressible Navier-Stokes equations using discrete boolean elements on a lattice. Rothman and Keller 
subsequently generalised the basic lattice gas method to allow simulation of immiscible fluids [ pof , and we 
have used their model as the starting point for our own work. 

Notwithstanding the simplifications engendered by invoking the lattice-gas paradigm, the simulation 
of the non-equilibrium behavior of ternary amphiphilic fluids in three dimensions is a highly demanding 
area of computational science; indeed, the results presented in this paper have been made possible only 
by the recent availability of sufficiently powerful parallel computing architectures, as well as sophisticated 
visualisation methods. 

The purpose of the present paper is to describe the implementation of our three-dimensional model, and 
to establish its validity. In particular, we show that our model can reproduce the well known features of 
amphiphilic fluid phenomenology. In Sections ^| and [}] we describe our model, emphasising in particular 
the differences between the 2D and 3D lattice-gas realisations, and briefly describe the computational re- 
quirements of the work. Section || outlines the basic structure of the algorithm, while Section || specifies 
the coupling constants which are used in our simulations. The results of the simulations are presented in 
Section |^. These simulations demonstrate the ability of our model to represent binary immiscible behaviour, 
binary water-surfactant self-assembly and ternary amphiphilic behaviour. Finally we close the paper with 
some conclusions in Section \A 



2 Amphiphilic Lattice-Gas Dynamics 

Our lattice-gas model of amphiphilic fluids consists of three different species of particles moving about on 
a Z?-dimensional lattice C in discrete time steps. The three species are the two immiscible fluids (oil and 
water) denoted by red (R) and blue (B) colours, respectively, and the amphiphile A. Lattice-gas particles 
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can have velocities Cj, where 1 < i < b, and b is the number of velocities per site. We shall measure discrete 
time in units of one lattice update, so that a particle emerging from a collision at site x and time t with 
velocity c, will stream to site x + Cj where it may undergo the next collision. The Ci must thus be integer 
multiples of the lattice vectors; it is also possible to have = for some i to allow for "rest particles" with 
zero speed. 

We let nf (x, t) £ {0, 1} denote the presence (1) or absence (0) of a particle of species a £ {R, B, A} with 
velocity Cj, at lattice site x £ C and time step t. The nf (x, t) are not all independent since an "exclusion 
rule" is enforced whereby there can be no more than one particle of a given velocity at a given lattice site 
at a given time. The collection of all nf (x, t) for 1 < % < b shall be denoted by n a (x, t). This is not to be 
confused with the total number of particles of a given colour, 

b 

«"(x,i)^<(x,i). (1) 
»=i 

Likewise, we shall sometimes need the colour charge associated with a given site, 

qi {^t)=nf (x,t)-nf (x,t), (2) 
as well as the total colour charge at a site, 

b 

q(K,t) = Y,Qi{^,t). (3) 

i=l 

Finally, the collection of all n Q (x, t) for a £ {R,B,A} will be called the population state of the site; it is 
denoted by 

n(x,t) £ N (4) 

where we have introduced the notation TV" for the (finite) set of all distinct population states. 

In addition to the specification of the particle populations at a site, the amphiphile particles have an ori- 
entation, denoted by <Ji (x, t). This orientation vector, which has fixed magnitude a, specifies the orientation 
of the director of the amphiphile particle at site x and time step t with velocity c^. (Of course, if there is no 
amphiphile particle with that site, time step and velocity then the value of er^ (x, t) there is not defined.) In 
our work, the values of the er^ (x, t) vectors may vary continuously on a (D — l)-sphere, denoted by S ^ 1 ; 

It is also possible to construct models with discrete set of values for the CTj, but we do not consider that possibility in this 
paper. 
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thus, they take their values on a circle for D = 2 and a sphere for D = 3. The collection of the b vectors 
<Ji (x, t) at a given site x and time step t is called the orientation state, 

b 

£(x,i) = (o-i(x,i),<r 2 (x,i),...<7 i ,(x,i)) 6 ®S D -\ (5) 

where <8> b denotes the 6-fold Cartesian product. This is not to be confused with the total director of a site 

fa 

<x (x, t)=J2 <n t) . (6) 

i=l 

We shall also find it useful to define the following scalar director field 

fa 

S(x,t)=J2c t -(T t (x,t). (7) 

i=l 

Together, the population state and the orientation state completely specify (in fact, as noted above, they 
ovcrspecify) the total state of the site, 

s=(n,£), (8) 

where we have omitted the site and time-step specification for brevity. 

The time evolution of the system is an alternation between a streaming or propagation step and a collision 
step. In the first of these, the particles move in the direction of their velocity vectors to new lattice sites. 
This is described mathematically by the replacements 

n?(x + Ci,t+l)«-n?(x,t) (9) 

o- i (x + c i ,i + l)<-<T i (x,t), (10) 

for all x G C, 1 < i < b and a G {i?, B, A}. That is, particles with velocity Cj simply move from point x to 
point x + Cj in one time step. 

In the collision step, the newly arrived particles interact, resulting in new momenta and directors. The 
collisional change in the state of a lattice site x is required to conserve the mass of each species present 

b 

p a (x,^]T<(x,t), (11) 

i 

as well as the /^-dimensional momentum vector 

fa 

p(x,^^r]Tc 4 <(x,i) (12) 



(where we have assumed for simplicity that the particles all carry unit mass). Thus, the set TV of population 
states at each site is partitioned into equivalence classes of population states having the same values of these 
conserved quantities. For a site with given masses p = (p R , p B , p A ) and momentum p, we denote the set of 
allowed population states by E(p, p) C TV. Since the conservation laws do not restrict the orientations of 
the directors, the set of allowed total states is then 

b 

£(p,p)=E(p,p)<g)S D - 1 . (13) 

Given a precollision total state s e E(p, p) with masses p and momentum p, the postcollision total state 
s' must belong to the set 

s'= (n',S') ef(p(«),p(«)). (14) 

Henceforth, primed quantities are understood always to refer to the postcollision state. In our model, s' 
is sampled from a probability density V(s'), sometimes equivalently denoted T^n'jS'), imposed upon this 
set. We assume that the characteristic time for collisional and orientational relaxation is sufficiently fast in 
comparison to that of the propagation that we can model this probability density as the Gibbsian equilibrium 
corresponding to a Hamiltonian function; that is 

V(s') = ^exp[-l3H(s')}, (15) 

where f3 is an inverse temperature, H(s') is the energy associated with collision outcome s', and Z is the 
equivalence-class partition function, 



z( P , p, /?) = I d£ ' cx p \-p h w\ , 

n'C W.( 11 nl J 



(16) 



n'£B(p,p) 

and we have defined the measure on the set of orientational states 

b 



J = (g) J dai. (17) 

i— 1 



In practice, we sample s' = (n', S ) by first sampling the postcollision population state n' from the reduced 
probability density 

P(n') = J dT,'V{s'). (18) 



We then sample the postcollision orientation state by sampling the b orientations crj from each of 




for 1 < i < b; these are, as we shall see, independent distributions, so the b samples may each be taken 
without regard for the other outcomes. This completes the collision process. 



3 Local Amphiphilic Lattice-Gas Hamiltonian 

It remains to specify the local Hamiltonian used in the collision outcome selection process. Such a Hamilto- 
nian has been derived and described in detail for the two-dimensional version of the model @, and we use 
the same one here Q It is 

H(s') = 3 ■ (aE + //P) + o> ■ (eE + CP) + J ■ (e£ + (P) + -v(x, t) 2 , (20) 

where we have introduced the colour flux vector of an outgoing state 

b 

J(x,t) =^c 4 q ? '(x,t), (21) 

i=l 

the dipolar flux tensor of an outgoing state 

b 

J{x,t) = ^2atr' i (x,t), (22) 



! = 1 



the colour field vector 



the dipolar field vector 



the colour field gradient tensor 



E(x,t) =J2 c il(^ + Ci,t), (23) 



P(x,i) = -^CiS(x + Ci,t), (24) 



£ (x,t) =^ Ci E(x + c 4 ,i), (25) 



i=i 



2 Actually, the forms given for the various fields in reference |g are somewhat more general than those given here. In this 
presentation we restrict ourselves to nearest neighbour interactions for simplicity. 
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the dipolar field gradient tensor 



b 



(26) 



i=l 



and the kinetic energy of the particles at a site. 




(27) 



where v is the average velocity of all particles at a site, the mass of the particles is taken as unity, and a, /i, 
e, C an d S are coupling constants. 

We note that the change in kinetic energy was not included by Rothman and Keller in their immiscible 
lattice gas model. We include it here for completeness, and set S = 1.0; although it will not affect the 
equilibrium properties of the model, it may impact the non-equilibrium properties. It should be further 
noted that the inverse temperature-like parameter (3 is not related in the conventional way to the kinetic 
energy. For a discussion of the intoduction of this parameter into lattice gases to reproduce critical behaviour 
we refer the reader to the original work by Chan and Liang j^]. 

Eqs. ((2^ - [26]) were derived by assuming that there is an interaction potential between colour charges, 
and that the directors are like "colour dipoles" in this context The terms containing a models the 
interaction of colour charges with surrounding colour charges as in the original Rothman-Keller model [ pof ; 
those containing [i model the interaction of colour charges with surrounding colour dipoles; those containing 
e model the interaction of colour dipoles with surrounding colour charges (alignment of surfactant molecules 
across oil- water interfaces); those containing Q model the interaction of colour dipoles with surrounding 
colour dipoles (interfacial bending energy or "stiffness"). 

Note that the field quantities depend only on the precollision state, whereas the flux quantities depend 
on the postcollision state. Thus, the fields may be computed once at every site, just after the propagation 
step. The fluxes, on the other hand, must be computed for each possible outgoing state. It follows that the 
Hamiltonian may be written in the form 



3 Note that this definition differs from that used in the reference [Q. The two definitions can be reconciled by a redefinition 
of the coupling constants. 



h 




(28) 
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where we have defined 

Ho(n') = J • (aE + fjiP) (29) 

and 

A, = a [eE + CP + (e£ + CP) ■ c *l ■ (30) 
The reduced probability density for n' is then 

P(n') = / dEfVtf = CXPh/ f o(n0] fl / d < CX P (~^f^ ■ °<) , (31) 

i— 1 

and this demands evaluation of integrals of the form 

J dcr cxp (-n/3A • er) . (32) 

In two dimensions (D = 2), we can adopt a polar coordinate 0, chosen so that 6 = is the direction of 
A, so that this integral becomes 

/■27T 

/ dO cxp (-npa\ A| cos (9) = 27r/ (n/3a| A|) , (33) 
Jo 

where Iq is the modified Bessel function of the first kind |2J. So, from from Eq. Jig), we get the reduced 
probability density 

p|il)g w,wwi ^ ih , Wii|| (34) 

1=1 

which must be evaluated numerically for every possible outgoing state. 

In three dimensions (D = 3), we adopt spherical coordinates with polar axis in the direction of A, so 
that the integral over the orientational degrees of freedom is of the form 

/■2-7T pTT 1 if n = 

/ d(j) d9 smO exp (-n/3cr|A| cos 9) = 4ir < (35) 

J0 Jo sinh^ojAjj :f„_i 

[ /3 CT |A| un-i, 

where is the colatitude and <fi is the azimuthal angle. The reduced probability density of Eq. ( |l8|) is then 



,. (47r) b exp[-/3H (n')] Afi , „Ju 



n, 



amh(J3a\Ai\) 
pa\Ai\ 



(36) 



i=l 

Once the reduced probability is sampled from P(n'), we turn our attention to the determination of the 
new orientation state. For D = 2, the reduced probability density for the angle ()[ is given by 

(a ,s cxp(-/3g|A, t |cos6>;) 

* m = ^WW) ■ (37) 

10 



For D = 3, the probability density for the colatitude 9[ and azimuthal angle </>J is then 

cxp (— /3<r [Aj| cos 0-) sin 



47T 



.ih(;3er|Aj|) 



/3a|A, 



(38) 



4 Algorithmic Considerations 

A computer implementation of our hydrodynamic lattice-gas model requires a choice of data representation 
for the population state n and the orientation state S. We consider these in turn. As noted above, the 
variables nf (x, t) arc not all independent because of the "exclusion rule" that only one particle of any type 
may have a given velocity i at a given site x and time step t. Thus, it is inefficient to store these variables 
directly. Rather, we note that for a given i, x and t there are precisely four possibilities: there can be a 
particle of type R, type B, type A, or nothing at all. These four possibilities can be encoded in two bits of 
information as follows: 



High Bit 


Low Bit 


Description 








Nothing 





1 


R particle 


1 





B particle 


1 


1 


A particle 



Thus, the population state of a given site can be represented by 2b bits of information. For D = 2, our 
current implementation uses a triangular lattice (coordination number 6) and one rest particle, so b = 7 
and 14 bits are needed to store the population state [Q. For D = 3, our current implementation uses a 
projected face-centered hypercubic (PFCHC) lattice (coordination number 24) and two rest particles, so 
b = 26 and 52 bits are needed to store the population state [0,^3|. The PFCHC lattice will be described in 
some detail later in this section. For now we note that in either case, the population state easily fits in a 
single integer variable; more precisely, for D = 2 it fits in a "short" integer of 16 bits, and for D = 3 it fits 
in a "double-precision" integer of 64 bits. 
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The orientation state requires much more storage because we have chosen to allow the orientation angles 
to be continuous ^. For D = 2 we must store the b polar angles each as an IEEE-format, single precision, 
floating-point variable of 32 bits. For D = 3 we must likewise store the two spherical angles for each velocity 
for a total of 2b floating-point variables. While it is true that these angular variables arc not defined unless 
the corresponding nf' = 1, our current implementation provides for storage for all angles at all sites and 
velocities, because the computational price of dynamically allocating and deallocating variables was not 
thought to be worth the savings in storage; for very low surfactant concentrations, this assumption may be 
invalid, and so the existing code might be improved. 

In a language that allows for user-defined types, such as Fortran 90, C++ or Java, it is useful to create 
a type for the total state of a site that includes both the population and orientation information. Given 
this data representation, the implementation of the propagation step is fairly obvious. The substitutions 
of Eqs. (Q) and (|l0|) are made throughout the lattice. In Fortran 90, the CSHIFT intrinsic accomplishes 
periodic shifts on arrays, and it is natural to use this to construct a subroutine that accepts an array of 
type 'total state' and performs the propagation procedure on this array as a side effect. The above-described 
representation for the population state is somewhat inconvenient in this regard, as the bit pairs corresponding 
to a particular velocity must be extracted from the integer variable before it is shifted in that direction. When 
the propagation step is completed, the new fields S*, E, P, £ and V are computed at each site using Eqs. (Q) 
and © - ©. 

Next, the possible outcomes for the population state are enumerated using a lookup-table procedure. A 
list of all possible outgoing states that has been sorted according to equivalence class E(p, p) is prccomputed 
and stored. This list is of length 2 14 = 16384 for the D = 2 case. A full list for the 52-bit population state 
representation in D = 3 would have length 2 52 and this is obviously much too large to store on any existing or 
contemplated computer; a method of shortening this list will be described below. For now, assume that such 
a list could be stored, and that two other lists of the same size could be maintained that accept the current 
population state n as a key, and return a pointer to the initial position and the number of elements of the 

4 As noted in an earlier footnote, this choice is not thought to be essential, and it is possible that much storage space may 
be saved by requiring the orientation angles to be discrete. 
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equivalence class E(p, p) in the table of sorted states. This makes it possible to enumerate the postcollision 
states n' that respect the required conservation laws. Note that the length of this list may vary from site to 
site. 

Next, each site loops over its set E(p, p) of allowed postcollision states and computes the outgoing colour 
flux vector J and the b auxiliary vectors Aj for each one of them, using Eqs. ( pl| ) and ( |30| ) respectively. 
These are then used to compute the reduced probability density P(n'), given by Eq. ( |34| ) for D = 2, or 
Eq. ( |36| ) for D = 3. Given these probabilities, a final population state n' is sampled. 

Once n' is known, the Aj are recalculated for that final state, and the final orientation angles are sampled 



from Eq. (37) for D = 2, or Eq. ( B8[ ) for D = 3. In the latter case, we note that TTi(8' i ,4>'i) is independent 
of <^ so this may simply be sampled uniformly in (0, 2tt). The colatitude 9[ is then found by equating its 
cumulative distribution function to a random number r uniformly distributed between and 1, and solving 
for 9'-. The result is 



-1 



In 



Ja\A t \ 

For D = 2, the sampling procedure is not so simple, because of the integral leading to the modified Bessel 
function. In this case we proceed numerically; for small values of the parameter /3cr|Aj|, we approximate 
the distribution by a Gaussian centered at its maximum, and for small values of (3a\Ai\ we employ rejection 
sampling [Q]. 

It remains to describe the PFCHC lattice and our method for making the size of the lookup table more 
manageable. The face-centered hypercubic (FCHC) lattice is a regular, self-dual lattice in four dimensions 
with coordination number 24. It can be described as all integer-valued tetrads k, I) such that i+ j + k + l 
is even. The motivation for using this lattice is that it is known to yield isotropic Navier-Stokes behaviour 
for a single-phase fluid . 

The lattice vectors are then the 24 neighbours of the site (0,0,0,0), and these can be partitioned into 
a subset of eight lattice vectors called Group 1, namely (±1,0,0, ±1) and (0, ±1,±1,0), a subset of eight 
lattice vectors called Group 2, namely (0, ±1, 0, ±1) and (±1, 0, ±1, 0), and a subset of eight lattice vectors 
called Group 3, namely (0, 0, ±1, ±1) and (±1, ±1, 0, 0). The virtue of this partition is that the sixteen lattice 
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vectors of any two groups lie on the corners of a regular four-dimensional hypercube, and the eight lattice 
vectors of the remaining group point to the face centers of that hypercube Q . 

The projection of the FCHC lattice to the three-dimensional PFCHC lattice can be accomplished by 
simply ignoring the fourth coordinate of the 24 vectors described above. The partition into three groups 
of eight vectors is still useful to maintain, as we shall see. One feature of this projection is that distinct 
vectors of the FCHC lattice can project to the same vector on the PFCHC lattice; for example, (1, 0, 0, 1) and 
(1,0,0, —1) both project to (1,0,0). We take these 24 three-vectors as the particle velocities in our D = 3 
model, and add two rest particles to them for a total of 26 particle velocities (b = 26). In our computer 
implementation, we append the (same) two rest particles to each of our three groups of eight lattice vectors, 
to obtain three groups of ten velocities (b — 10) each. The idea is then to allow collisions within each group 
of ten velocities separately, updating the state of the rest particles in all three groups whenever they change, 
thereby letting the rest particles provide mass and momentum transfer between the three groups. The three 
sets of velocities are summarized in the following table: 



Group 


Lattice 


Components 


Group 


Lattice 


Components 


Group 


Lattice 


Components 




Vector 


X, 


y, 


z 




Vector 


X, 


y, 


z 




Vector 


X, 


y, z 




Cl 


1, 


o. 







Cl 


o, 


i, 







Cl 


o, 


0, 1 




c 2 


1, 


o. 







c 2 


o, 


i, 







c 2 


o, 


0, 1 




c 3 


-1, 


0, 







c 3 


o, 


-i, 







c 3 


o, 


0, -1 




c 4 


-1, 


0, 







c 4 


o, 


-i, 







c 4 


o, 


0, -1 


1 


c 5 


0, 


i. 


1 


2 


c 5 


1, 


0. 


1 


3 


c 5 


1, 


i, o 




c 6 


0, 


i. 


-1 




c 6 


-1. 


0. 


1 




c 6 


1, 


-1, 




C7 


o, 


-i, 


1 




C7 


1, 


0, 


-1 




C7 


-1, 


1, 




c 8 


o, 


-i, 


-1 




c 8 


-1, 


0, 


-1 




c 8 


-1, 


-i, o 




c 9 


o, 


0. 







c 9 


o, 


0. 







eg 


o, 


0, 




ClO 


o, 


0. 







ClO 


o, 


o, 







ClO 


o, 


0, 
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Since two bits of information are required to represent the population state for each velocity, a total of 20 
bits will specify the state within each group. This results in tables of length 2 20 = 104 8 5 76. Since the results 
can be stored as single-precision (32-bit, or four-byte) integers, the tables each require a total of 4 Mbytes 
of storage. Since there are three tables, as described above, a grand total of 12 Mbytes are devoted to table 
storage. This amount of storage is not a significant problem on modern multiprocessor supercomputers. 

Once again, the use of user-defined types can simplify the implementation of the above decomposition. 
The population state can be made a user-defined type, constructed from three integer variables. Subroutines 
for propagation, table lookup, and other collision-related operations can then be overloaded to accommodate 
the new type. 

Despite the advantage gained by using the hydrodynamic lattice gas method described above, the simu- 
lation of large three dimensional amphiphilic systems remains extremely computationally demanding. The 
algorithm described above has been implemented in Fortran90, and parallelised using the MPI (Message 
Passing Interface) library. Simulations for the parameter search described in Section ^| were performed on a 
512 processsor Cray T3D at the Edinburgh Parallel Computing Centre (EPCC). The code was then ported 
to SGI Origin 2000 machines at Schlumberger Cambridge Research, the National Computational Science 
Alliance (NCSA) in Illinois, and at the Oxford Supercomputing Centre. The Cray T3D and T3E systems 
are Massively Parallel Processor (MPP) machines, whereas the SGI O2000 machines have cache-coherent 
Non-Uniform Memory Architecture (ccNUMA). These two parallel architectures are very different, and the 
ease of moving from one to the other illustrates the portability of modern parallelised codes. A performance 
improvement of two orders of magnitude was obtained on going from the T3D to the O2000 machines, and 
all the simulations described in Section || were performed on the latter platforms. Baseline performance for 
a 64 3 system running on 8 processors is 3.3 timesteps per minute for a binary oil-water system. Performance 
scales superlinearly out to 64 processors for 64 3 and 128 3 systems. Computational complexity increases as 
L as one increases the linear system size L. Currently the largest attainable system size is 128 3 due to 
memory limitations of the O2000 machines currently available to us. However, 256 3 and 512 3 systems are 
attainable within the limits of the O2000 architecture. A port to the Cray T3E (MPP architecture) has also 
been performed, and a comparison on a processor-by-processor basis gives a performance of 30-50% of the 
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O2000, with linear as opposed to superlinear scalability. However, the larger number of processors available 
on a typical T3E system more than compensates for this relative fall off-in one case we had access to a 
1200-node T3E for a period of one month. A full description of the computational aspects of this model will 
be provided elsewhere. 

In addition to the demanding nature of the simulations themselves, visualising the results produced can 
be as-and in some aspects more-computationally demanding. In particular, the generation of geometrical 
information required to plot 3D isosurfaces of individual species concentration requires RAM resources of at 
least 1 Gbyte. Although visualisation is today sometimes still regarded as a subsidiary activity to numerical 
simulation, we have found it absolutely vital to check that the code was working, to distinguish different 
morphoplogics and to gain intuition about the very complex dynamics of these systems. Visualisation of the 
largest and most complex systems attainable by simulation using our hydrodynamic lattice gas demands use 
of the most advanced graphics engines currently available. 



5 Definition of Coupling Constants 

The Hamiltonian used to determine the collision outcomes has been specified in Eq. (p0[). We now describe 
the choice of the coupling comstants a, (i, e, (. In addition to the terms in Eq. ( pp[ ) there is an additional 
kinetic energy term, so that the Hamiltonian becomes: 

H(s') = J • (aE + fiP) + cr' ■ (eE + (P) + J : (e£ + (V) + ^ |v(x, t)\ 2 , (40) 

where v is the average velocity of all particles at a site, and the mass of the particles is set to unity 
(5 = 1.0). Our model reduces to the Rothman-Kcllcr model for binary immiscible fluids in the limit a — > oo. 
The three remaining parameters control the surfactant interactions. These were chosen by performing an 
extensive parameter search using binary water-surfactant systems, and measuring the structure factor for 
these systems to look for signs of self-assembly. In particular, we sought strong structure- factor peaks 
indicative of spherical or wormlike micelles of a characteristic size were sought. These simulations were 
performed with a surfactant: water ratio of 1:2, well above the critical concentration for the formation of 
micelles. The physical contributions of each term in Eq. (|2(]), and therefore the effect of varying each 
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parameter, is described below: 

• e controls the interaction of outgoing dipolcs with the surrounding colour field. In ternary phases this 
term will send surfactant to oil-water interfaces. In binary water-surfactant phases this interaction will 
tend to favour flat interfaces between the phases; 

• (j, controls the interaction of outgoing coloured particles with the surrounding dipolar field. This term 
will favour the bending of dipoles around a central colour charge, and will be important in creating 
micellar phases; 

• £ controls the surfactant-surfactant interaction. For positive £ this produces an attractive, ferro- 
magnetic interaction between the dipoles. This term is of limited importance in the formation of 
self- assembled phases, and was set to 0.5 for the simulations described below. 

The key to locating the micellar phases is to find the correct balance between e and (i. Strongly peaked 
structure functions were obtained for the following values of e and ji: 

0.75 < // < 2.0 
0.25 < e < 2.0 

In order to maximise the desired behaviour of sending surfactant to oil-water interfaces while retaining the 
necessary micellar binary phases we chose a canonical set of coupling constants which are kept constant 
throughout all the following simulations, except in the short vesicle study described in Section |6.5| . The 
values of these constants are: 

a = 1.0, e = 2.0, \i = 0.75, C = 0.5 
The temperature-like parameter, /?, was held constant at 1.0 for all of the ensuing simulations. 

6 Simulations 

The equilibrium properties and phase diagrams of a wide variety of real amphiphilic fluids are well known |L4| . 
These phase diagrams have many features unique to three dimensional systems. To esablish the general 
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validity of our model it is essential that we can reproduce the complex phase behaviour observed in real 
amphiphiles. In this section we describe some of this phenomenology, and present results of the simulations 
designed to reproduce it. 

6.1 Oil- Water System 

The spinodal decomposition of immiscible fluids has been extensively studied in two dimensions, and rather 
less in three ||,^]. After a quench into the two phase region an initially homogeneous mixture separates 
into relatively pure single-phase domains. 





Figure 1: Binary phase separation in an oil-water system. Timcstcps (clockwise from top left) 200, 400, 600, 
1000. Red isosurface shows interface between oil and water phases. System size is 64 3 . 



With no surfactant particles present in the system, the only term in the local site Hamiltonian, Eq. (EOj), 
that contributes numerically to the collision process is aAH cc . With parameters a = 1.0, (3 = 1.0 we 
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performed simulations starting from an initial configuration in which the lattice vectors at each site are 
populated randomly with oil or water particles with equal probability (a so-called critical quench). This 
initial condition corresponds to a homogenised 1:1 oil water mixture. The reduced density (i.e. the proportion 
of lattice vectors at each site that contain a particle) is 0.5. In this parameter regime, the model exhibits 
phase separation with positive surface tension, as is evident from Fig. |], which illustrates the nonequilibrium 
behaviour of the immiscible lattice gas. If left to run for a large enough time the system would eventually 
reach the completely separated state of two distinct layers of fluid. To make a detailed comparison between 
the immiscible oil- water fluid behaviour shown here and later simulations in which we introduce surfactant to 
the system, we make use of spherically averaged structure functions. We first calculate the three-dimensional 
structure factor of the colour charge, s(k, t), 



s(k,t) = 




where k = (2tt/L) (pi + qj + rk), p,q,r = 1,2, q(x) is the total colour charge at a site, q av is the 

average value of the colour charge, L is the length of the system and N = L 3 is the number of lattice sites 
in the simulation box. 

We actually compute the spherically averaged structure factor, given by 

S(k,t) = ^^, (41) 

where k = 2nn/L, n = 0, 1, 2, L, and the sum is over a spherical shell defined by (n — i) < < 
(n+ i). The resolution of S(k, t) depends on k c , the Nyquist cutoff frequency associated with the lattice (for 
a real-space sampling interval of A the cut-off frequency is 1/(2A)); above this value of the frequency there 
is only spurious information due to aliasing. In our case, k c = (27r/L)n c , where we have chosen n c to be the 
maximum possible value, which is half the lattice size. Fig. |^ shows the temporal evolution of S(k, t) for 
the case of two immiscible fluids. As time increases, the peak of S(k, t) shifts to progressively smaller wave 
numbers and its peak height increases. This behaviour is characteristic of domain coarsening and serves as 
a useful comparison for the surfactant-containing systems described below. 
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Figure 2: Structure factor for timesteps 200, 400, 600, 800, 1000 in a binary phase-separating system . 
System size is 128 3 . The structure factor measurements are averaged over ten independent simulations. 
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6.2 Binary amphiphilic fluid phases: from monomers to sponges. 

In two dimensions the behaviour of the binary water-amphiphile lattice gas fluid is characterised by a critical 
micelle concentration (CMC) below which the fluid is a solution of monomeric amphiphiles, and above which 
the monomers form circular micelles. Further increase in amphiphile concentration in 2D simply results 
in more micelles, the micelles retaining their characteristic size. The situation in three dimensions is more 
complex. A CMC is still present for the formation of spherical micelles, but at higher concentrations new 
structures appear. As the concentration of surfactant increases the number of spherical micelles increases, 
until a second critical concentration is reached, beyond which wormlikc micelles are the preferred struc- 
ture. When the concentration of surfactant is high enough that both surfactant and water phases percolate 
throughout the system a sponge phase results. We identify this sponge phase with the L 3 phase described 
by Gompper and Schick fL4f . These concentrations have been determined in our model for /? = 1 as: 

0.006 < Pspherical < 0.012 
0.08 < Pwormlike < 0.25 

0.25 < p Ls 

where p is the reduced density; that is, the fraction of lattice sites occupied by surfactant particles. The 
description of these different regimes as distinct phases, with 'critical' concentrations separating them, may 
be somewhat misleading. There is a large degree of phase coexistence. Individual monomers may join or 
leave a micelle in the spherical micelle regime, and the kinetics of simple micelle formation can be modelled 
theoretically on the basis of a Becker-Doring theory |8|,[l^]. The critical concentration for the formation of 
spherical micelles is well defined in our model, with no micelles seen below this concentration. The formation 
of more complex structures, however, appears to take place by more general (Smoluchowski-type) agregation 
processes; wormlike micelles are formed from the coalesence of spherical micelles. Such behaviour reflects 
the highly dynamic nature of our model. 

We performed simulations designed to access the monomeric, spherical micelle, wormlike micelle and 
sponge phases. All started with random initial conditions and the coupling constants defined in Section |^. 
To determine the stability of these structures, we calculated spherically-averaged structure functions of the 
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Figure 3: Spherical micelles in water at surfactant concentration 0.008, timestep 1000. Left figure shows a 
32 3 system. Arrows represent the average surfactant direction per site, overlaid on an isosurface of surfactant 
concentration. Right figure shows a 64 3 system. The left half of this figure displays arrows of the average 
surfactant direction at a site, while the right half shows an isosurface of surfactant concentration. 



Monomers 

■ — — Spherical micelles 
- — - Wormlike micelles 
L 3 Sponge phase 



400 600 
Time (timesteps) 



Figure 4: Average value of the magnitude of the wavevector A: as a function of time for varying surfactant 
concentrations on a 64 3 lattice. The monomeric surfactant solution (solid line) is at reduced surfactant 
concentration of 0.001; spherical micelles at reduced surfactant concentration of 0.01; wormlike micelles at 
reduced surfactant concentration of 0.1; and the L3 sponge phase at reduced surfactant concentration of 0.3 
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surfactant density. In Fig. [|, we plot the temporal evolution of the characteristic wave number, 

t k t t )) = ^k=o kS ( k ^) 

the inverse of which is a measure of the average domain size. We see that in the low amphiphile concentration 
case the characteristic size remains consistent with a random configuration of solubilised monomers. For the 
cases where visualisation establishes the presence of more complex structures there is fast initial growth of 
the surfactant domains which soon levels off to some constant size; see Figs || and ||. 

6.3 Ternary phases: lamellae 

We first investigate the stability of a lamellar structure, which is composed of alternating layers of oil-rich 
and water-rich phases separated by a monolayer of surfactant molecules. We look at the system with and 
without surfactant present in order for a comparison to be made. It is not clear what influence surfactant 
will have on such a lamellar structure in three dimensions. We arc interested in the ability of surfactant to 
stabilise large areas of interface by lowering the interfacial tension. We set up the initial configuration of 
the system, with layers of oil and water four sites wide, all sites having a reduced density of 0.5. It is clear 
that if our model is exhibiting the correct behaviour, then we would expect there to be a critical density of 
surfactant required at the oil-water interfaces in order for the layered structure to be stable. Consequently, 
we set up a simulation with a single site wide layer of surfactant at each of the oil- water interfaces. Snapshots 
from these simulations are shown in Figs. ^ and Fig. |^; the former is the pure oil- water case with coupling 
coefficient a = 1.0 and inverse temperature /3 = 1.0, while the latter has surfactant monolayers present. 

The fluctuations present in the model enable oil and water particles from the initially separated layers 
to move; they thus come under the influence of the colour field gradients produced by other layers of the 
same fluid type. For these parameter choices there is an inherent tendency for the oil and water to act as 
immiscible fluids and phase separate (c/., Fig. |l|), precisely as is shown in Fig. [j]. Fig. || shows the case 
where surfactant monolayers coat the oil-water interfaces. Here we see that the initial periodic structure is 
stabilised despite the presence of large amounts of oil-water interface. 
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Figure 5: Wormlike micelles in water at a reduced surfactant concentration of 0.1 for a 32 3 system. The 
same structures are present in larger systems, but we display this snapshot for clarity. The green isosurface 
shows the surfactant concentration at a level of 5 particles per lattice site. 
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Figure 7: Binary oil-water lamellae at timesteps t = (left) and t = 500 (right) for a 48 3 system. The red 
and blue colourings on the slice planes indicate oil and water respectively, while the blue isosurface shows 
the interface between oil and water. For clarity the isosurface is shown for only half of the system. The 
system size is 48 3 




Figure 8: Ternary lamellae at t = (left), and t = 1000 (right). The isosurface shows the oil- water interface. 
System size is 40 3 . 
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6.4 Ternary phases: oil-in-water (water-in-oil) and bicontinuous microemul- 





Figure 9: Oil-in-water droplet phase arising from a random initial condition. Oil (red, left) and surfactant 
(green, right) isosurfaces. 



The ternary behaviour of amphiphilic systems shows the same increase in complexity between two and 
three dimensions that we have seen in the binary water-surfactant systems. One simplifying feature of 
our model is the symmetry between oil and water in the interactions producing the phase behaviour. We 
discuss here the two most basic of the many ternary phases, the oil-in-water droplet and the bicontinuous 
microemulsion phases. 

In the oil-in-water droplet phase, oil exists as the minority phase in a continuous background of water 
and surfactant. If sufficient surfactant is present in the system the spinodal decomposition of a quenched 
uniform mixture of the three components will be arrested, and the system will reach an equilibrium state 
consisting of droplets of oil in water with surfactant stabilising the oil- water interfaces. This phase may 
be regarded as a perturbation of a binary spherical micelle phase, the micelles now being swollen with oil. 
If one increases the oil concentration, a regime is reached where both the oil and water domains percolate 
throughout the system. Such a phase is the ternary extension of the L3 sponge phase in the binary system. 

In both these cases, the complete separation of the oil and water phases, which is the equilibrium state 
for a binary immiscible fluid, is prevented by the presence of the surfactant. In order to quantify this result, 
and to compare with the previous binary case, we calculate the mean k value of the spherically averaged 



structure factor of the colour charge, and plot the result in Fig. 10 
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Figure 10: Structure factor mean k value for binary phase separation, ternary coalescing droplet phase, 
swollen wormlike micellar phase, and bicontinuous microcmulsion phases. 

In order to reproduce the oil-in-water droplet phase, we set up two simulations with a random initial 
configuration and a reduced density of oil of 0.05. The concentration of surfactant was varied between the 
two simulations, the reduced surfactant density being 0.01 and 0.21 respectively. The total reduced density 
for both simulations was 0.5; to maintain consistency between our various simulations, we continue to use 
the numerical values for the coupling constants defined in Section [s]. 

Fig. H displays a snapshot of the simulation with a reduced surfactant density of 0.01 at timestep 1000, 
displaying the oil (red) and surfactant (green) isosurfaces. It is clear that the surfactant is adsorbed at 
the interface. However, Fig. |l0| shows the time evolution of the average k value of the spherically averaged 
structure factor. It is clear that although the phase separation is slowed by the presence of surfactant, the 
average domain size in the system is increasing with time. For the case with reduced surfactant density 0.21, 
one observes complete arrest of phase separation, with wormlike oil domains remaining at a constant size as 
shown by the dashed line in Fig. [HJ We refer to this morphology as a swollen wormlike micellar phase. A 
snapshot of this simulation taken at timestep 4000 is shown in Fig. [H]. 
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Figure 11: Bicontinuous microcmulsion phase at timestep 4000. Left: water isosurfacc; middle: surfac- 
tant isosurface; right: Oil isosurfacc. The system size is 64 3 . Isosurfaces display concentrations of water, 
surfactant and oil, respectively, at a level of five particles per site. 

To obtain the bicontinuous sponge regime within the model's phase diagram, we simply increase the 
relative amount of oil present in the system. Hence this simulation has a random initial mixture with a 
reduced density 0.5 and a 0.83 : 1.0 : 0.83 oil-to-surfactant-to- water ratio. Fig. |l2| displays the results of this 
simulation at timestep 4000, displaying the percolating oil (red) and water (blue) isosurfaces, as well as the 
surfactant (green) isosurfacc. 




Figure 12: Stable oil- in- water droplet phase at timestep 4000. Red and green isosurfaces show oil and 
surfactant concentrations at a level of five particles per site. The system size is 64 3 . 

These results show the ability of our model to correctly reproduce the effects of surfactant on the phase 
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separation dynamics of binary immiscible fluids. A quantitative study of the domain growth behaviour for 
both binary and ternary systems will be described in a future paper. 



6.5 Vesicles 

Membrane theory treatments of amphiphilic systems predict that for large binary water surfactant systems a 



bilayer can reduce its energy by closing its boundary, creating a large structure known as a vesicle |14|. Such 
structures can exist because the energy required to bend the bilayer into a sphere is less than that required 
to maintain the bilayer edge. We do not expect vesicles to self-assemble in our current simulations, in part 
because such structures are usually metastable in real life and require the input of energy, for example via 
shearing or sonication. We can, however, construct such a structure as an initial condition and study its 
stability, or lack thereof, using our model. The initial condition we have used is shown in Fig. |l3|. The 
concentration of surfactant is initially zero except within a spherical shell, 5 sites wide and 32 sites in radius, 
where the reduced density of surfactant is 1.0. All other sites contain water particles with a reduced density 
of 0.5. The total system size is 128 3 . 

For the following simulations, we have not used the canonical choice of coupling constants described in 
Section]^. Rather, in order to study the mechanisms by which these vesicles disintegrate, we hsve performed 
three separate simulations, with the following parameter choices: 

a = 1.0, e = 0.01, /i = 0.01, C = 1-0 
a = 1.0, e = 0.01, /i = 1.00, C = 0.01 
a = 1.0, e = 1.00, n = 0.01, C = 0.01 

These three simulations should separate the effects of the three surfactant interactions on the vesicles. We 
also wish to characterise the timescale over which the vesicles are stable. In order to do this, we performed 
one simulation with oil particles replacing the surfactant particles in the initial condition, and with a — 0.0. 
This essentially labels the particles starting in the vesicle region, and allows us to track their subsequent 
motion. Not suprisingly, this simulation shows that the initial vesicle configuration is destroyed in the 
absence of interactions in less than ten timcstcps. The initial condition and state of the system at timestep 
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200 for the above three simulations is shown in Fig. O. These simulations show that the mechanisms driving 
the breakup of the vesicles are phase separation and micellisation. The vesicles are also unstable against 
fluctuations in concentration which produce holes in the structure. Once a hole has been produced it is 
observed to grow until the vesicle is destroyed. 
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Figure 13: Vesicle disintegration. Clockwise from top left: Initial condition. Simulation with e = 0.01, 
/i = 1.00, C = 0.01 at timcstcp 200. Simulation with e = 1.00, n = 0.01, £ = 0.01 at timcstcp 200. 
Simulation with e = 0.01, fj, = 0.01, ( = 1.0 at timcstcp 200. All snapshots display an isosurface showing 
reduced surfactant density at a level of 5 particles per site for half the vesicle. 
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7 Conclusions 

We have extended our hydro-dynamic lattice-gas model for the dynamics of binary and ternary amphiphilic 
fluids from two to three dimensions. We have shown that our model exhibits the correct 3D phenomenology 
using a combination of visual and analytic techniques. Experimentally observed self-assembling structures 
form in our simulations in a consistent manner when the relative concentrations of the three components are 
varied. Binary immiscible, binary amphiphilic and ternary amphiphilic behaviour are all captured using a 
single set of coupling constants. We have also shown that studies of vesicle metastability arc possible using 
this model with different choices of the coupling constants. Work is currently in progress on a wide range of 
amphiphilic fluid systems, including measurements of viscosity and surface tension, and the study of growth 
laws in amphiphilic self-assembly processes. Studies of amphiphilic fluid flow in porous media, which have 
previously been performed in 2D |Q , are now underway using the 3D model [|| . Indeed there is a rich seam 
of problems related to amphiphilic fluids which may be mined using this model. Our work confirms the 
suitability of lattice gas automata for the modelling and simulation of such complex fluid problems in both 
two and three dimensions. 
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